function [averageBetas maskedBetas] = extractROImask(subjID, ADAPT_t_thresh, FFAcontrast, subjectsDir)
%
%   _______________________________________________
%   by Marcelo G Mattar (12/03/2012)
%   mattar@sas.upenn.edu


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% CHECK INPUTS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if nargin < 4
    subjectsDir = '/Users/marcelomattar/Data/AdaptID/Subjects/';
end

if nargin < 3
    FFAcontrast = 2; %1 for Faces-Objects, 2 for Faces-Scenes, 3 for Main Effect
end

if nargin < 2
    ADAPT_t_thresh = 1.65;
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% BASIC DEFINITIONS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


% Define some basic directories and files
ROIsdir = [subjectsDir subjID '/ROIs/'];
shortTermFaceDir = [subjectsDir '/' subjID '/ShortTermFace/mbs5f50_1adaptCov/'];
indexOf_1back_adaptation_within_feat = 3;


% Retrieve the analysis directories
FaceDirectories_feat = struct2cell(dir([shortTermFaceDir 'Face*.feat']));
FaceDirectories_feat = FaceDirectories_feat(1,:)';
numFaceRuns = length(FaceDirectories_feat);


averageBetas = cell(numFaceRuns+2,2);
averageBetas{1,1} = subjID;

averageBetas{2,1} = 'Run #1';
averageBetas{3,1} = 'Run #2';
averageBetas{4,1} = 'Run #3';
if numFaceRuns == 4
    averageBetas{5,1} = 'Run #4';
end
averageBetas{numFaceRuns+2,1} = 'Average';

averageBetas{1,2} = '1-back Adaptation';
averageBetas{1,3} = 'Standard Deviation';
averageBetas{1,4} = 'Mask size';

maskedBetas = cell(numFaceRuns,1);
averageBetasMatrix = zeros(numFaceRuns,1);
stdBetasMatrix = zeros(numFaceRuns,1);
maskSizes = zeros(numFaceRuns,1);

% Pre-allocate some important variables
feat_cope = cell(numFaceRuns,1);
leaveOneOut_tstat = cell(numFaceRuns,1);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% LOAD FFA MASK
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if exist([ROIsdir 'lFFA' num2str(FFAcontrast) 'mask.nii.gz'],'file')
    lFFAmask = loadniigz([ROIsdir 'lFFA' num2str(FFAcontrast) 'mask.nii.gz']);
else
    possibleContrasts = [1 2];
    lFFAmask = loadniigz([ROIsdir 'lFFA' num2str(possibleContrasts(~(FFAcontrast == [1 2]))) 'mask.nii.gz']);
end

if exist([ROIsdir 'rFFA' num2str(FFAcontrast) 'mask.nii.gz'],'file')
    rFFAmask = loadniigz([ROIsdir 'rFFA' num2str(FFAcontrast) 'mask.nii.gz']);
else
    possibleContrasts = [1 2];
    rFFAmask = loadniigz([ROIsdir 'rFFA' num2str(possibleContrasts(~(FFAcontrast == [1 2]))) 'mask.nii.gz']);
end

FFAmask = lFFAmask + rFFAmask;
FFAmask = logical(FFAmask);






%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% LOAD ADAPT DATA
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Extract Adaptation copes, varcopes and tstats for each of the Face runs
for runIndx = 1:numFaceRuns
    feat_cope{runIndx,1} = loadniigz([shortTermFaceDir FaceDirectories_feat{runIndx} '/reg_standard/stats/cope' num2str(indexOf_1back_adaptation_within_feat) '.nii.gz']);
end

for leftout = 1:numFaceRuns
    leaveOneOut_tstat{leftout,1} = loadniigz([shortTermFaceDir 'Face_LeaveOneOut/Face_LeaveOut' num2str(leftout) '.gfeat/cope' num2str(indexOf_1back_adaptation_within_feat) '.feat/stats/tstat1.nii.gz']);
    thisMask = and((leaveOneOut_tstat{leftout,1} > ADAPT_t_thresh), FFAmask);
    maskSizes(leftout,1) = sum(thisMask(:));
    maskedBetas{leftout,1} = feat_cope{leftout,1}(thisMask(:));
    averageBetasMatrix(leftout,1) = mean(maskedBetas{leftout,1});
    stdBetasMatrix(leftout,1) = std(maskedBetas{leftout,1});
end






%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% SAVE RESULTS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


for i=1:numFaceRuns
    averageBetas{i+1,2} = averageBetasMatrix(i,1);
    averageBetas{i+1,3} = stdBetasMatrix(i,1);
    averageBetas{i+1,4} = maskSizes(i,1);
end

averageBetas{numFaceRuns+2,2} = mean(averageBetasMatrix);
averageBetas{numFaceRuns+2,3} = std(averageBetasMatrix);
averageBetas{numFaceRuns+2,4} = ' ';
